/*
 * Author:
 *	Pierre Lindenbaum PhD.
 * Contact:
 *	plindenbaum@yahoo.fr
 * WWW:
 *	http://plindenbaum.blogspot.com
 * Motivation:
 *	Generates a Postscript Manhattan Plot 
 * Date:
 *	2011
 */
#include <cstdlib>
#include <vector>
#include <map>
#include <cerrno>
#include <string>
#include <cstring>
#include <stdexcept>
#include <climits>
#include <cmath>
#include <cfloat>
#include <cstdio>
#include <iostream>
#include <zlib.h>
#include <sstream>
#include <algorithm>
#include <cassert>
#include <stdint.h>
using namespace std;

#define THROW(a) do{ostringstream _os;\
	_os << __FILE__ << ":"<< __LINE__ << ":"<< a << endl;\
	throw runtime_error(_os.str());\
	}while(0);

typedef int32_t pos_t;
typedef float pixel_t;
typedef double value_t;
#define TEST_LT(a) if(this->a!=cp.a) return this->a < cp.a

struct Date
	{
	int year;
	int month;
	int day;
	bool operator < (const Date& cp) const
		{
		TEST_LT(year);
		TEST_LT(month);
		return day < cp.day;
		}
	};

struct AbstractData
    {
	value_t y;
    };

struct XYData:public AbstractData
	{
	value_t x;
	bool operator < (const XYData& cp)
		{
		return x < cp.x;
		}
	};

struct TemporalData:public AbstractData
	{
	Date x;
	bool operator < (const TemporalData& cp)
		{
		return x < cp.x;
		}
	};

enum DataType
	{
	TEMPORAL,XY
	};

class DataFactory
	{
	public:
		virtual size_t sizeOfData()const=0;
		virtual bool parse(const char* src,void* dest)const=0;
		virtual bool lt(const void* a,const void* b)const=0;
		virtual string toString(const void* c)const=0;
		virtual pixel_t convert(const void* min,const void* max, const void* value,pixel_t width)=0;
		virtual void sort(void* data, size_t count) const;
	};

class XYDataFactory:public DataFactory
	{
	private:
			value_t cast(const void* p)
				{
				return *((value_t*)p);
				}

	public:
			virtual size_t sizeOfData()const
				{
				return sizeof(value_t);
				}
			virtual bool parse(const char* src,void* dest)const
				{
				char* p2;
				value_t v=(value_t)strtod(src,&p2);
				if(*p2!=0) return false;
				memcpy((void)dest,&v,sizeof(value_t));
				}
			virtual bool lt(const void* a,const void* b)const
				{
				return cast(a) < cast(b);
				}
			virtual string toString(const void* c)const
				{
				ostringstream os;
				os << cast(c);
				return os.str();
				}

			virtual pixel_t convert(const void* min,const void* max, const void* value,pixel_t width)
				{
				return (cast(value)-cast(min))/(cast(max)-cast(min))*width;
				}
	private:
			static int compare (const void * a, const void * b)
				{
				value_t d=cast(a)-cast(b);
				return d==0?0:d<0?-1:1;
				}
	public:
			virtual void sort(void* data, size_t count) const
				{
				qsort(data,count,sizeOfData(),_compare);
				}
	};


class Serie
	{
	public:
		DataFactory* factory;
		string name;
		void* data;
		void* min;
		void* max;
		size_t buffer;
		size_t size;
		Serie(DataFactory* factory):factory(factory)
			{
			data=NULL;
			buffer=0;
			size=0;
			min=NULL;
			max=NUL;
			}
		~Serie()
			{
			free(data);
			free(min);
			free(max);
			}
		void add(const void* newval)
			{
			if(size+1>=buffer)
				{
				buffer=(buffer==0?BUFSIZ:buffer*2);
				data=safeRealloc(data,buffer*factory->sizeOfData());
				std::memcpy(
						&((char*)data)[factory->sizeOfData()*size],
						newval,
						factory->sizeOfData()
						);

				if(min==NULL)
					{
					min=safeMalloc(factory->sizeOfData());
					memcpy(min,newval,factory->sizeOfData());
					}
				else if(factory->lt(newval,min))
					{
					memcpy(min,newval,factory->sizeOfData());
					}

				if(max==NULL)
					{
					max=safeMalloc(factory->sizeOfData());
					memcpy(max,newval,factory->sizeOfData());
					}
				else if(factory->lt(max,newval))
					{
					memcpy(max,newval,factory->sizeOfData());
					}

				size++;
				}
			}

		pixel_t convert(const void* value,pixel_t length)
			{
			return factory->convert(min,max,value,length);
			}
		void sort()
			{
			return factory->sort(data,size);
			}

	};

class Plot
    {
    private:
		bool readline(gzFile in,std::string& line)
			{
			line.clear();
			int c;
			if(gzeof(in)) return false;
			while((c=gzgetc(in))!=EOF && c!='\n')
				{
				line+=(char)c;
				}
			return true;
			}
    public:
		DataFactory* vFactory;
		DataFactory* hFactory;
		Serie* vSerie;
		vector<Serie*> hSeries;
		char delim;

		std::ostream* output;
		pixel_t x_axis_width;
		pixel_t y_axis_height;
		pixel_t margin_left;
		pixel_t margin_right;
		pixel_t margin_top;
		pixel_t margin_bottom;

		Plot():delim('\t'),
			chrom2info(),
			output(&cout)
			{
			margin_left=100;
			margin_bottom=100;
			margin_top=100;
			margin_right=100;
			x_axis_width=841.89 -(margin_left+margin_right);
			y_axis_height=595.28-(margin_top+margin_bottom);
			}

		~Plot()
			{
			for(vector<Serie*>::iterator r=hSeries.begin();
				r!=hSeries.end();
				++r)
				{
				delete (*r);
				}
			}

		void split(const string& line,vector<string>& tokens)
			{
			tokens.clear();

			}

	void readData(gzFile in)
	    {
		char *hData=malloc(hFactory->sizeOfData());
		char *vData=malloc(hFactory->sizeOfData());
	    string line;
	    vector<string> header;
	    if(first_line_is_header)
	    	{
	    	if(!readline(in,line))
	    		{
	    		cerr << "Cannot read 1st line / header" << endl;
	    		exit(EXIT_FAILURE);
	    		}

	    	split(line,header);
	    	vSerie=new Serie(vFactory);
	    	vSerie->name.assign(header[0]);
	    	for(size_t i=1;i< header.size();++i)
	    		{
	    		Serie* serie=new Serie(hFactory);
	    		serie->name.assign(header[i]);
	    		hSeries.push_back(serie);
	    		}
	    	}
	    size_t nLine(0);
	    vector<string> tokens;
	    while(readline(in,line))
			{
	    	nLine++;
			if(line.empty()) continue;
			split(line,tokens);
			if(!first_line_is_header)
				{
				vSerie=new Serie(vFactory);
				vSerie->name.assign("$1");
				for(size_t i=1;i< tokens.size();++i)
					{
					ostringstream os;
					os << "$"<<(i+1);
					Serie* serie=new Serie(hFactory);
					serie->name.assign(os.str());
					hSeries.push_back(serie);
					}
				}
			else if(tokens.size()!=header.size())
				{
				cerr << "Expected "<< header.size() << " columns but got "
						<< tokens.size()
						<< endl;
				exit(EXIT_FAILURE);
				}
			if(!vFactory->parse(tokens[0].c_str(),vData))
				{
				cerr << "Error line "<< nLine
					<< " Cannot parse \"" << tokens[0] << "\""<< endl;
				exit(EXIT_FAILURE);
				}
			vSerie->add(vData);

			for(size_t i=1;i< tokens.size();++i)
				{
				if(!vFactory->parse(tokens[i].c_str(),hData))
					{
					cerr << "Error line "<< nLine
						<< " Cannot parse \"" << tokens[0] << "\""<< endl;
					exit(EXIT_FAILURE);
					}
				hSerie[i-1]->add(vData);
				}
			}
	    free(hData);
	    free(vData);
	    }
	pixel_t pageWidth()
	    {
	    return (margin_left+margin_right+x_axis_width);
	    }
	pixel_t pageHeight()
	    {
	    return (margin_top+margin_bottom+y_axis_height);
	    }
	void print()
	    {


	    for(chrom2info_t::iterator r=chrom2info.begin();
		    r!=chrom2info.end();
		    ++r)
		{
		ChromInfo* c=r->second;
		if(c->chromStart==c->chromEnd) c->chromEnd++;
		double length5=1+(c->length()/100.0)*5.0;
		c->chromStart=max(0,(pos_t)(c->chromStart-length5));
		c->chromEnd=(pos_t)(c->chromEnd+length5);

		std::sort(c->data.begin(),c->data.end());

		size_of_genome+=c->length();
		assert(size_of_genome>0);
		}
	    pixel_t x=margin_left;
	    for(chrom2info_t::iterator r=chrom2info.begin();
	   		    r!=chrom2info.end();
	   		    ++r)
		{
		ChromInfo* c=r->second;
		c->x=x;
		c->width=(c->length()/(double)size_of_genome)*x_axis_width;
		assert(c->length()>0);
		assert(size_of_genome>0);
		assert(x_axis_width>0);
		assert(c->width>0);
		x+=c->width;
		}

	    ostream& out=(*output);
	    time_t rawtime;
	    time ( &rawtime );
	    out <<
		    "%!PS\n"
		    "%%Creator: Pierre Lindenbaum PhD plindenbaum@yahoo.fr http://plindenbaum.blogspot.com\n"
		    "%%Title: " << __FILE__ << "\n"
		    "%%CreationDate: "<< ::ctime (&rawtime)<<
		    "%%BoundingBox: 0 0 "
			<<(int)pageWidth()<< " "
			<<(int)pageHeight()<<"\n"
		    "%%Pages! 1\n"
		    "/marginBottom "<< margin_bottom << " def\n"
		    "/marginLeft "<< margin_left << " def\n"
		    "/yAxis "<< y_axis_height << " def\n"
		    "/getChromName { 0 get } bind def\n"
		    "/getChromStart { 1 get } bind def\n"
		    "/getChromEnd { 2 get } bind def\n"
		    "/getChromX { 3 get } bind def\n"
		    "/getChromWidth { 4 get } bind def\n"
		     "/getChromLength {1 dict begin /chrom exch def chrom getChromEnd  chrom getChromStart sub end} bind def\n"
		    "/getChromDataIndex { 5 get } bind def\n"
		    "/getChromDataCount { 6 get } bind def\n"
		    "/convertPos2X {2 dict begin /pos exch def /chrom exch def  "
		    "pos chrom getChromStart sub chrom getChromLength div chrom getChromWidth mul chrom getChromX add\n"
		    "end } bind def\n"
		    "/minValue "<< minValue<<" def\n"
		    "/maxValue "<< maxValue<<" def\n"
		    "/toString { 20 string cvs} bind def\n"
		    "/Courier findfont 14 scalefont setfont\n"
		    "/convertValue2Y {\n"
		    " minValue sub maxValue minValue sub div yAxis mul marginBottom add"
		    "} bind def\n"
		    "/dd { 2 dict begin\n"
			"  /y exch def\n"
			"  /x exch def\n"
			" newpath "
			"   x y 5 add moveto\n"
			"   x 5 add y lineto\n"
			"   x y 5 sub lineto\n"
			"   x 5 sub y lineto\n"
			" closepath 0 0  1 setrgbcolor fill "
			"   end } bind def\n"
		    "/box { 4 dict begin\n"
		    "  /height exch def\n"
		    "  /width exch def\n"
		    "  /y exch def\n"
		    "  /x exch def\n"
		    "   x y moveto\n"
		    "   width 0 rlineto\n"
		    "   0 height rlineto\n"
		    "   width -1 mul 0 rlineto\n"
		    "   0 height -1 mul rlineto\n"
		    "   end } bind def\n"
		    ;
	    int index_data=-1;
	    for(chrom2info_t::iterator r=chrom2info.begin();
		    r!=chrom2info.end();
		    ++r)
		    {
		    index_data++;
		    ChromInfo* chromInfo=r->second;
		    
		    out << " "<< (index_data%2==0?1:0.7)<< " setgray newpath "<< chromInfo->x << " "<< margin_bottom << " "<< chromInfo->width << " " << y_axis_height <<  " closepath box fill ";
		    
		    out << " 0.4 setgray newpath "<< chromInfo->x << " "<< margin_bottom << " "<< chromInfo->width << " " << y_axis_height <<  " closepath box stroke ";
		    }
		    
		    
		for(chrom2info_t::iterator r=chrom2info.begin();
		    r!=chrom2info.end();
		    ++r)
		    {
		    ChromInfo* chromInfo=r->second;
		    for(vector<Data>::iterator r2=r->second->data.begin();
			    r2!=r->second->data.end();
			    ++r2)
			    {
			    out << chromInfo->convertToX(&*(r2))  << " " << convertToY(&(*r2)) << " dd\n";
			    index_data++;
			    }
		    if(chromInfo->width>200)
		    	{
		    for(int i=0;i<=10.0;++i)
	   		{
	   		double x=chromInfo->x + ((chromInfo->width)/10.0)*i;
	   		out << " 0 setgray newpath "
	   		    << " " << x << " "  << (margin_bottom) << " moveto "
	   		    << " 0 -10 rlineto closepath stroke " 
	   		    << " " << x << " "  << (margin_bottom-20) << " moveto "
	   		    << " -90 rotate (" << (int)(chromInfo->chromStart+((chromInfo->length()/10.0)*i)) << ") show 90 rotate"
	   		    << endl
	   		    ;
	   		}
	   		}
	   	    
	   	    out 	<< " 0 setgray "
	   	    		<< " " << (chromInfo->x + (chromInfo->width/2.0))
	   	    		<< " "  << (margin_bottom+10+y_axis_height) << " moveto "
	   		    	<< "  90 rotate (" << chromInfo->chrom << ") show -90 rotate "
	   		    	<< endl;
		    }

	
	   for(int i=0;i<=10;++i)
	   	{
	   	out << " 0 setgray newpath ";
	   	out << " " << margin_left << " " << (margin_bottom+(i/10.0)*y_axis_height) << " moveto -15 0 rlineto ";
	   	out << " closepath stroke ";
	   	out << " 1 " << (margin_bottom+(i/10.0)*y_axis_height) << " moveto";
	   	out << " (" << (minValue+ (i/10.0)*(maxValue-minValue)) << ") show "; 
	   	}
	  
	   
	    out << "0 setgray\n0 0 "<< pageWidth()<<" " << pageHeight() << " box\n"
		"0 setgray\n"
		"stroke\n"
		;
	    out << "showpage\n";
	    out.flush();
	    }

    };

int main(int argc,char** argv)
    {
    Manhattan app;
    int optind=1;
    while(optind < argc)
   		{
   		if(std::strcmp(argv[optind],"-h")==0)
   			{
   			cerr << argv[0] << "Pierre Lindenbaum PHD. 2011.\n";
   			cerr << "Compilation: "<<__DATE__<<"  at "<< __TIME__<<".\n";
   			exit(EXIT_FAILURE);
   			}
   		else if(std::strcmp(argv[optind],"-h")==0)
   			{
   			cerr << argv[0] << "Pierre Lindenbaum PHD. 2011.\n";
   			cerr << "Compilation: "<<__DATE__<<"  at "<< __TIME__<<".\n";
   			exit(EXIT_FAILURE);
   			}
   		else if(argv[optind][0]=='-')
   			{
   			fprintf(stderr,"unknown option '%s'\n",argv[optind]);
   			exit(EXIT_FAILURE);
   			}
   		else
   			{
   			break;
   			}
   		++optind;
                }
    if(optind==argc)
	    {
	    gzFile in=gzdopen(fileno(stdin),"r");
	    if(in==NULL)
		{
		cerr << "Cannot open stdin" << endl;
		return EXIT_FAILURE;
		}
	    app.readData(in);
	    }
    else
	    {
	    while(optind< argc)
		{
		char* filename=argv[optind++];
		gzFile in=gzopen(filename,"r");
		if(in==NULL)
		    {
		    cerr << "Cannot open "<< filename << " " << strerror(errno) << endl;
		    return EXIT_FAILURE;
		    }
		app.readData(in);
		gzclose(in);
		}
	    }
    app.print();
    return EXIT_SUCCESS;
    }
